Rayleigh-Ritz Calculation of Effective Potential Far From Equilibrium 
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We demonstrate the utility of a Rayleigh-Ritz scheme re- 
cently proposed to compute the nonequilibrium effective po- 
tential nonperturbatively in a strong noise regime far from 
equilibrium. A simple Kramers model of an ionic conductor 
is used to illustrate the efficiency of the method. 

PACS. Numbers: 02.50.-r, 05.40. +j 



We have recently proposed to use in nonequilibrium 
statistical mechanics a novel variational principle associ- 
ated to the so-called effective action functional JlJ . This 
quantity is well-known in quantum field theory, where 
the concept has it roots in the early work of Heisenberg 
& Euler H and Schwinger Q in QED. In nonequilib- 
rium statistical mechanics, the first such action principle 
seems to have been Onsager's 1931 "principle of least dis- 
sipation" ||, which applies to systems subject to thermal 
or molecular noise, governed by a fluctuation-dissipation 
relation. A formulation of the least-dissipation principle 
by an action functional on histories was later developed 
by Onsager and Machlup 1|. In a weak- noise limit the 
field-theoretic effective action and the Onsager-Machlup 
action coincide, as discussed some time ago by Graham 
0. (See also ||). However, the weak-noise limit is of 
rather restricted applicability. In particular, it is useless 
to deal with problems in which there is no small param- 
eter and in which strong fluctuations dominate the phe- 
nomena on a wide range of length-scales. These include, 
for example, high Reynolds-number fluid and plasma tur- 
bulence, spinodal decomposition, and surface growth by 
random deposition. In strong- noise systems, efficient cal- 
culational tools remain to be discovered. One powerful 
nonperturbative scheme which was developed in quan- 
tum field theory is a Rayleigh-Ritz method based upon a 
constrained variation of the quantum energy-expectation 
functional p-pj|. This method has been recently ex- 
tended by us to nonequilibrium statistical dynamics with 
non-Hermitian evolution operator . 

It is our purpose here to demonstrate the computa- 
tional utility of that method in a simple concrete model, 
the Kramers model of an ionic conductor fll2]. The 



model consists of a unit mass, charged particle in a 1- 
dimensional cosine potential V(x) — Vocosir, damped 
by linear friction with coefficient 7, and driven both by 
an electric force F = qE and by thermal noise of tem- 
perature T — 70/fcs. The dynamics of individual real- 
izations is given by the nonlinear Langevin equation 



jx + V Q smx = F + T(t), 



with zero-mean Gaussian noise: 



{T(t)T(t')) = 2jeS(t-t'). 



(1) 



(2) 



Here x runs over [0, 2ir] with periodic b.c. Equivalently, 
the evolution of the single-time distribution P t of posi- 
tion x and velocity v = x is given by the Fokker-Planck 
equation 
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(vPt) + — [{jv + V cosx 

7 Q^-)P t ] = LP t . 



(3) 



The important point for our discussion here is that all 
parameter values will be chosen order one, i.e. Vq — 
7 = F = 6 = 1. In this case, there is neither a small 
nor a large parameter on which to base a perturbation 
expansion or asymptotic evaluation. In particular, F = 1 
implies that the system is far from thermal equilibrium, 
well beyond the regime of linear response of steady-state 
current j(F) = (v) p to the applied field F. 

In the statistical steady state, rather than the full 
effective-action, it is more useful to consider its time- 
extensive limit, the so-called effective potential. For any 
random variable z(t) in an ergodic system, this quan- 
tity measures the "cost" for fluctuations to occur in the 
empirical time-averages: 



dt z(t). 



(4) 



That is, it measures the probability for deviations to oc- 
cur in a given realization at large T between zt and the 
ensemble-mean z. Quantitatively, the relation between 
this probability and the effective potential function V is 



1 



that, for any possible single-time value z of the variable 
z(f), the probability for the average over a time-interval 
of duration T to take on the value z has magnitude 



Prob ({z r « z}) - exp (-T ■ V[z}) 



(5) 



in the limit of large T. The potential function V[z] is non- 
negative and convex, with minimum value =0 occuring 
only at the ensemble-mean, V[z] = 0. The fluctuation 
formula Eq.(|^) is a refinement of the standard ergodic 
hypothesis. It states not only that time-averages zt will 
converge to the ensemble-average z almost surely in the 
limit of large T, but also it gives a precise quantitative 
estimate on the probability of the deviations. This refine- 
ment holds whenever some condition of finite exponential 
moments is satisfied: see The effective potential is 

analogous to the thermodynamic entropy of an equilib- 
rium system (more precisely, to its negative) and Eq.(|^) 
for the probability of statistical deviations is analogous 
to the Einstein-Boltzmann formula for the fluctuations 
away from equilibrium. This probabilistic interpretation 
of the effective potential was pointed out in field-theory 
by Jona-Lasinio [jO). Although any variable might be 
considered, we shall be interested here in the electric cur- 
rent j and its effective potential V\f\. 

The method we use to calculate V is the Rayleigh- 
Ritz variational scheme proposed in As discussed 

in more detail there, the effective potential may be char- 
acterized very generally by means of a constrained vari- 
ation of the "energy functional" 



(6) 

defined in terms of the Liouville operator L of the 
nonequilibrium statistical dynamics. V[z] was shown to 
be the extremal value of 7i\%> , \E ,i ] under arbitrary vari- 
ation of "trial states" ty R > L subject to the constraints of 
constant overlap 



(* L ,# fl ) = 1 



and constant means 



(7) 



(8) 



The operator Z consists of multiplication by variable z. 
By means of Lagrange multipliers A, h associated to the 
constraints, this variational characterization of the effec- 
tive potential can be expressed in spectral terms as 



V[z] = z ■ h A[h], 



(9) 



where A[h] is the principal eigenvalue of the "perturbed" 
Liouville operator = L + h • Z. Observe that this 
spectral characterization is exactly analogous to the char- 
acterization of the free-energy in equilibrium statistical 
mechanics as the principal eigenvalue of the "transfer- 
matrix" in lattice systems (13]. In that case, the en- 
tropy is then obtained by a Legendre transform of the 



free-energy, entirely analogous to the relation of A[h] 
and V[z] in Eq.(^|) by a Legendre transform. As the 
starting point for an approximate Rayleigh-Ritz calcula- 
tion of V, any probability distribution function (PDF) 
ansatz for the stationary distribution of the nonequilib- 
rium dynamics may be employed. The "right trial state" 
^(x) just consists of a guess for the stationary PDF, 
p(x;c), depending upon some arbitrary parameters c. 
The "left trial state" \t (x) is a (linear superposition of 
a) set of selected moment- functions (x) . The variation 
within such a framework leads in general to a "nonlinear 
eigenvalue problem" to determine the approximate value 
Ajt[h] of the leading eigenvalue within the ansatz. See 



The simplest systematic approximation procedure in 
the Kramers model is to use the Gaussian trial weight 



w(x, v) 



1 
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(10) 



along with the corresponding truncated orthogonal poly- 
nomial expansions 



JY 



'(*,„) =«,(*, to-x; x (ii) 

n=0p=—P v v 



and 
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Here He n (v) are a set of Hermite polynomials, with con- 
ventions as in |fl5|| . The same expansion was used by 
Risken & Vollmer |l(J to calculate the full stationary 
measure. Orthogonal polynomial expansions have the 
simplifying feature that they lead to a linear problem 
for the approximate eigenvalue A#[/i] and eigenvectors 
c # j c # (with shorthand for N, P) , as 



and 



with 



X( L #. fl )"P<"'P' C #,« , P' = ^#[%#,np> ( 13 ) 
n' ,p' 
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(15) 



for n = 0,...,N,p = -P,...,P and (Jy# th )nv,n'p' = 
CL#,h)n'p'.n P ' tne Hermitian conjugate. L#,h is the finite 
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matrix approximation to the full "perturbed" operator 
Lh = L + hV. The approximate value of the effective 
potential is then obtained from 



V # [j]=j # [h}-h-X # [h), 



(16) 



in which A# [h] is the eigenvalue branch passing through 
at h = and j#[h] = V^^{h)) is the associated 

"perturbed" current. The latter may be determined from 
the Hellmann-Feynman theorem jl7), as = A^[/i]. 

It can also be calculated directly from the approximate 
right and left eigenvectors, as 



E 

It. p 



#,n-l,p) 



l( c #,n+l,p)* c #,n,p] ■ (17) 



This avoids the evaluation of derivatives and is more com- 
putationally efficient in such a small system, where the 
determination of the eigenvalues and eigenvectors is easy 
to accomplish. These calculations have been carried out 
by us numerically, and the results for the approximate 
effective potential V# [j] are graphed in Figure 1 for vari- 
ous choices of N and P. As can be seen from that figure, 
convergence is already obtained in the range j = 0.2 — 1.4 
to 1% accuracy for N = 4, P = 5. 

The effective potential can also be obtained by di- 
rect numerical solution (DNS) of the Langevin dynam- 
ics. The most obvious procedure would be to gather 
long time series of v(t) in the steady-state, to assem- 
ble histograms of probabilities of empirical time averages 
T>t = Jq dt v(t) over time- intervals of duration T, and 
then to calculate V[j] via the inverse to Eq.(||). However, 
we have found that the approximate potential from this 
method converges very poorly. Instead, the most efficient 
method is to determine approximate values of A# [h] and 
j#[h] from 



and 



X#[h] = ^ln(exp 



3#[h] 



h I dt v(t) 
o 



(18) 
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where (•) denotes average over some number R of real- 
izations and # now stands for R, T and the numerical 
time step At. Substituting these into Eq.(|l6|), a value 
V#\j] is obtained which for large R, T and small At con- 
verges to the true potential. We have also performed 
this direct calculation of V[j] via a 2nd-order stochas- 
tic Runge-Kutta integration of the Langevin Eq.(|) for 
R = 32000, T = 1000 and At = 0.01. Increasing T or 
decreasing At did not change the results, for which the 



largest errors are statistical, based upon the finite num- 
ber of realizations R. 

In Figure 2 we compare also the potentials V[j] them- 
selves for the above DNS and the N = 10, P = 8 
Rayleigh-Ritz, over the range j = 0.718 — 0.733. Over 
this range the error in the Rayleigh-Ritz value is at most 
one part in 10 4 . The relatively large error bars on Vdns 
are due to the cancellation of linear terms in h between 
\#[h] and j#[h] ■ h in Eq.(|l6|). Because of this cancel- 
lation, getting an accuracy even to one part in 10 2 for 
V# required an accuracy to one part in 10 4 for A#[/i] 
and j#[h]. To achieve an accuracy of Vdns compara- 
ble to Rayleigh-Ritz would require a reduction in sta- 
tistical error in DNS values of A# and j# from 1 part 
in 10 4 to 1 part in 10 8 . With this error estimated as 
0(i? -1 / 2 ), it would require a computation time longer by 
a factor of 10 8 . The present calculation required ~ 10 12 
floating-point operations (flops) and, on the machine we 
employed, took about 10 hrs, the bulk of which was de- 
voted to generating random numbers. Thus, accuracy to 
one part in 10 4 for V# from DNS would require about 
10 5 years of computation! By contrast, the Rayleigh- 
Ritz computation with iV = 10, P = 8 accurate to 0.01% 
required ~ 10 6 flops, performed in about 0.1 seconds. 
In general, the Rayleigh-Ritz evaluation of V#\j] in the 
Kramers model requires for each value of h a calcula- 
tion of the leading eigenvalue \#[h] and its associated 
eigenvector for the tridiagonal matrix Eq. (|l5|) , which is 
an 0(N ■ P) operation, and the calculation of current 
j#[h] via the summation formula Eq.(p^), which is also 
0{N ■ P). The number of flops required is - 10007V • P. 
The superiority of the Rayleigh-Ritz method compared 
with DNS is clear. Even more favorable is the fact that 
the Rayleigh-Ritz method gives a result as in Figure 1 
over a range 100 times larger, to 1% accuracy, with only 
TV = 4, P — 5. In this range errors in DNS due to finite T 
and R are very large because the /i-weighted ensembles 
in Eqs.(|l8|),(|l9|) have a greater contribution from rare 
events as h increases. Outside of a narrow range of cur- 
rent j near the mean value the determination of effective 
potential by DNS is practically impossible. 

In this work we have focused upon the Rayleigh-Ritz 
determination of the effective potential via a convergent 
scheme, the expansion in orthogonal polynomials. How- 
ever, as discussed in detail in O any physically-inspired, 
nonlinear ansatz for the PDF or any "surrogate" ran- 
dom variables chosen to model the system variables may 
be used as well. This is important for systems with 
many degrees-of-freedom, such as high Reynolds num- 
ber turbulence or large-scale dynamics of multiphase flu- 
ids, in which convergent schemes are presently not re- 
motely feasible. In work in progress [^8[ , we apply our 
variational methods to practical modeling of such sys- 
tems, in particular hydrodynamic turbulence. Analogous 
time-dependent Rayleigh-Ritz methods are available to 
calculate the full effective action functional on histories: 
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see |2|. There is some general similarity of our methods 
to the Hartee-Fock variational approximation applied by 
Crisanti & Marconi to the calculation of effective 
action in phase segregation dynamics. Both are intrinsi- 
cally nonperturbative and capable of systematic improve- 
ment. In addition, the Rayleigh-Ritz scheme illustrated 
here has a very great flexibility to incorporate intuitive 
guesses into an analytical calculation: any PDF ansatz 
for the system variables whatsoever may be used as the 
basis for a calculation by our variational method. 
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FIGURE CAPTIONS 

Figure 1.) Approximate effective potential V[j] for 
N=4, P=5 (□), N=6, P=6 (+) and N=10, P=8 (o). 

Figure 2.) Potentials V[j] for DNS (with errorbars) 
and the N = 10, P = 8 Rayleigh-Ritz solid line). 
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